Event-by-event evaluation of the prompt fission neutron spectrum from ^^^Pu(n, /) 
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Earlier studies of ^^®Pu(n, /) have been extended to incident neutron energies up to 20 MeV 
within the framework of the event-by-event fission model FREYA, into which we have incorporated 
multichance fission and pre-equilibrium neutron emission. The main parameters controlling prompt 
fission neutron evaporation have been identified and the prompt fission neutron spectrum has been 
analyzed by fitting those parameters to the average neutron multiplicity V from ENDF-B/VIl.O, 
including the energy-energy correlations in ^{E) obtained by fitting to the experimental V data 
used in the ENDF-B/VII.O evaluation. We present our results, discuss relevant tests of this new 
evaluation, and describe possible further improvements. 



I. INTRODUCTION 

Nuclear fission forms a central topic in nuclear physics, 
presenting many interesting issues for both experimen- 
tal and theoretical research, and it has numerous prac- 
tical applications as well, including energy production 
and security. Nevertheless, a quantitative theory of fis- 
sion is not yet available. While there has been consid- 
erable progress in the last few years, both in liquid-drop 
model-type calculations [H, 0] and in microscopic treat- 
ments [J-l^l, these treatments primarily address "cold" 
fission, induced by thermal neutrons, and cannot yet de- 
scribe "hot fission" , induced by more energetic neutrons. 
In order to perform new evaluations of observables impor- 
tant for applications over the full relevant energy range, 
it is therefore necessary to rely on a considerable degree 
of phenomenological modeling. 

One of the most important quantities for applications 
is the prompt fission neutron spectrum (PFNS). As dis- 
cussed earlier ^g], the experimental spectral data them- 
selves are neither sufficiently accurate nor of sufficiently 
consistent quality to allow an improved PFNS evaluation. 
However, by combining measured information about the 
nuclear fragment yields and energies with the very pre- 
cise evaluations of neutron multiplicities, it is possible to 
constrain the neutron spectrum rather tightly without 
having to rely on the spectral data themselves. 

Our approach employs the fission model FREYA (Fis- 
sion Reaction Event Yield Algorithm) which incorporates 
the relevant physics and contains a few key parameters 
that are determined by cornparison to pertinent data 
through statistical analysis [1, 01- It simulates the en- 
tire fission process and produces a large sample of com- 
plete fission events with full kinematic information on the 
emerging fission products and the emitted neutrons and 
photons. It thus incorporates the pre-fission emission of 
neutrons from the fissile compound nucleus as well as 
the sequential neutron evaporation from the fission frag- 
ments. FREYA provides a means of using readily-measured 
observables to improve our understanding of the fission 
process and it is, therefore, a potentially powerful tool 



for bridging the gap between current microscopic mod- 
els and important fission observables and for improving 
estimates of the fission characteristics important for ap- 
plications. 

In the following, we briefiy describe the employed ver- 
sion of FREYA and the fitting procedure used to obtain 
our extended evaluation of the ^^^Pu(n, /) prompt fis- 
sion neutron spectrum. We then compare our results 
to the ENDF-B/VII.O [1] evaluation of the PFNS and 
some benchmark criticality tests. Finally, we discuss the 
energy and model dependence of several relevant obser- 
vaables. 



II. GENERATION OF FISSION EVENTS 

We have adapted the recently developed fission model 
FREYA for the present purpose of calculating the neu- 
tron spectrum in terms of a set of well-defined model 
parameters. We describe its main physics ingredients be- 
low, with an emphasis on the new features added for the 
present study, particularly multichance fission and pre- 
equilibrium emission. Being a simulation model, FREYA 
follows the temporal sequence of individual fission events 
from the initial agitated fissionable nucleus, 240p^* 
the present case, through possible pre-fission emissions 
to a split into two excited fragments and their subse- 
quent sequential emission of neutrons and photons. The 
description below is similarly organized. 



A. Pre-fission neutron emission 

At low incident neutron energies, below a few MeV, 
the neutron is absorbed into the target nucleus resulting 
in an equilibrated compound nucleus which may have a 
variety of fates. Most frequently, in the present case, it 
will fission directly. But, since the compound nucleus was 
formed by neutron absorption, it is energetically possible 
for it to re-emit a neutron. In that circumstance, the 
daughter nucleus cannot fission and will de-excite by se- 
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FIG. 1: (Color online) The probability for first- (black circles), 
second- (red squares), third- (green diamonds), and fourth- 
(blue triangles) chance fission as a function of incident neutron 
energy. The solid curves show the GNASH results used in the 
ENDF-B/VII.O evaluation [l], while the dashed and dot-dot- 
dashed curves with closed and open symbols are the FREYA 
results discussed in the text. 



quential photon emission. FREYA generally discards such 
events because it is designed to provide fission events (but 
their frequency is noted). Neutron evaporation from a fis- 
sionable compound nucleus can be treated in the same 
manner as neutron evaporation from fission fragments, 
as will be described later (Sect. IiTD ip . In principle, it is 
also possible that the compound nucleus will start by ra- 
diating a photon but the likelihood for this is very small 
and is ignored. 



Bf is the height of the fission barrier (the corresponding 
quantity for neutron emission is the neutron separation 
energy S„). 

Neutron evaporation is possible whenever the excita- 
tion energy of the compound nucleus exceeds the neu- 
tron separation energy, E* > Sn- (Since it costs energy 
to remove a neutron from the nucleus, 5„ is positive.) 
The excitation of the evaporation daughter nucleus is 
= E* — Sn — E where E is the kinetic energy of the rel- 
ative motion between the emitted neutron and the daugh- 
ter nucleus. If this quantity exceeds the fission barrier in 
the daughter nucleus, then second-chance fission is pos- 
sible. (We use the Hill- Wheeler expression for the trans- 
mission probability, Pj — 1/[1 + exp(27r(i3/ — E*j)/huj)\ 
with tujj—1 MeV, so there is an exponentially small prob- 
ability for sub-barrier fission.) The procedure described 
above is then applied to the daughter nucleus, thus mak- 
ing further pre-fission neutron emission possible. Thus as 
the incident neutron energy is raised, the emission of an 
ever increasing number of pre-fission neutrons becomes 
possible and the associated fission events may be classi- 
fied as first-chance fission (when there are no pre-fission 
neutrons emitted), second-chance fission (when one neu- 
tron is emitted prior to fission), and so on. 

Figure [T] shows the probabilities for n"^-chance fission 
obtained with FREYA for incident neutron energies up to 
20 MeV, using two alternate assumptions about the level 
density parameterization. Also shown are the GNASH re- 
suhs used in the ENDF-B/VII.O evaluation Q. The two 
calculations give rather similar results but, because these 
probabilities are not easy to measure experimentally, it is 
not possible to ascertain the accuracy of the calculations. 



Multichance fission 



2. Pre-equilibrium neutron emission 



As the energy of the incident neutron is raised, neu- 
tron evaporation from the produced compound nucleus 
competes ever more favorably with direct (first-chance) 
fission. The associated probability is given by the ratio 
of the fission and evaporation widths Ff (£'*) and Tn{E*) 
for which we use the transition-state estimate H, 



TniE*) _ 2gnHnCr J^^iXn- x)pn{x)dx 
T{{E*) ~ TTh? 



J^' pi{x)dx 



(1) 



where Qs — 2 is the spin degeneracy of the neutron, is 
its reduced mass, and a — ttR^ — tttqA^I'^. Furthermore, 
Pn {x) is the level density in the evaporation daughter nu- 
cleus at the excitation energy x, whose maximum value 
is given by A"„ = Qn = E* — Sn, where Qn is the Q value 
for neutron emission and S'„ is the neutron separation 
energy. Similarly, p[{x) is the level density of the tran- 
sition configuration for the fissioning nucleus, i.e. when 
its shape is that associated with the top of the fission 
barrier; the excitation x is measured relative to that bar- 
rier top, so its maximum value is Xf = E* — B{, where 
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FIG. 2: (Color online) Upper panel: The probability for pre- 
equilibrium neutron emission as a function of the incident neu- 
tron energy. Lower panel: The corresponding average multi- 
plicity of neutrons emitted prior to fission calculated without 
(dashed) and with (solid) the pre-equilibrium processes. 
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At higher incident neutron energies, there is a grow- 
ing chance that complete equihbrium is not estabhshed 
before the first neutron is emitted. Under such cir- 
cumstances the calculation of statistical neutron evap- 
oration must be replaced by a suitable non-equilibrium 
treatment. A variety of models have been developed 
for this process (for example Ref. Tq*] which combines 
pre-equilibrium emission with the Madland-Nix model 
of the prompt fission neutron spectrum) and we em- 
ploy a practical application of the two-component exciton 
model [HI (described in detail in Ref. 13]). It repre- 
sents the evolution of the nuclear reaction in terms of 
time-dependent populations of ever more complex many- 
particle-many-hole states. 

A given many-exciton state consists of p^(Tr) neutron 
(proton) particle excitons and h^^^-^ neutron (proton) 
hole excitons. The total number of neutron (proton) ex- 
citons in the state is n^(^Tr) — Pv{Tr) + h-i.^-^y The incident 
neutron provides the initial state consisting of a single ex- 
citon, namely a neutron particle excitation: = 1 and 
Pn = = = 0. In the course of time, the number 
of excitons present may change due to hard collisions or 
charge exchange, as governed by the residual two-body 
interaction. We ignore the unlikely accidental processes 
that reduce the number of excitons, so the state grows 
ever more complex. 

The temporal development of the associated probabil- 
ity distribution P{pu, /ii/jPtt, ^tt) is described by a master 
equation that accounts for the transitions between differ- 
ent exciton states. The pre-equilibrium neutron emission 
spectrum is then given by 
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where ctcn is the compound nuclear cross section (usu- 
ally obtained from an optical model calculation), W is 
the rate for emitting a neutron with energy E from 
the exciton state {pTr,hTr,Pi,,hi,), r is the lifetime of 
this state, and P(p^, /i,y) is the (time-averaged) 
probability for the system to survive the previous stages 
and arrive at the specified exciton state. In the two- 
component model, contributions to the survival proba- 
bility from both particle creation and charge exchange 
need to be accounted for. The survival probability for 
the exciton state (ptt, ^7r,Piy, ^i^) can be obtained from 
a recursion relation starting from the initial condition 
P{p,y = l,hi, = OjPtt = 0, /itt = 0) = 1 and setting P = 
for terms with negative exciton number. As in Ref. [l3| . 
particle emission is assumed to occur only from states 
with at least three excitons, n,^ + n^^ > 3. We consider 
excitons up to p'^'^^ = p'^'^^ — 6. 

The emission rate, W{pT^,hT^,Pi,,h^,Ek), is 
largely governed by the particle-hole state density, 
^{p-m h.„,Pv,h^, E*). For a neutron ejectile of energy E 
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FIG. 3: (Color online) The contributions to the pre- 
equilibrium neutron spectrum from exciton states with the 
indicated values of (pir, h^,pu,h^), obtained at En — 14 MeV. 
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Ul{p-n, K,Py - 1, K, E* - E - Sn) 

u{p^,K,Py,hy,E*) 



(3) 



where (Tfc^inv is the inverse reaction cross section (calcu- 
lated within the optical model framework) and E* is the 
total excitation energy of the system. 

The calculated probability for pre-equilibrium neutron 
emission is shown in the upper panel of Fig. [2] as a func- 
tion of the incident neutron energy £"„, while Fig.[3]shows 
the the pre-equilibrium neutron spectrum obtained at 
En = 14 MeV. After being practically negligible below 
a few MeV, the probability for pre-equilibrium emission 
grows approximately linearly to about 24% at 20 MeV. 
A careful inspection of the calculated energy spectrum 
shows that neutrons emitted from states with larger ex- 
citon number approach the statistical emission expected 
from a compound nucleus, thus ensuring our treatment 
has included sufficient complexity to exhaust the pre- 
equilibrium mechanism. Because of the (desired) insen- 
sitivity to the maximum specified exciton number, the 
probability shown in Fig.[2]is not indicative of the impor- 
tance of the pre-equilibrium processes. Their quantita- 
tive significance is better seen by comparing the neutron 
spectrum obtained with and without the pre-equilibrium 
treatment, as shown in the lower panel of Fig. [2j 

The reaction cross sections used in Eqs. © and ^ 
define the overall magnitude of the cross sections for 
the pre-equilibrium processes. The highest accuracy re- 
sults are best obtained from coupled-channels calcula- 
tions with an appropriately-determined optical potential. 
However, since FREYA principally deals with probabilities, 
the relative fraction of pre-equilibrium neutrons may be 
computed with sufficient accuracy employing a spherical 
optical potential to calculate the relevant cross sections. 
Consequently, the compund-nucleus cross sections and 
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inverse cross sections were computed using the optical- 
model program ECIS06 and the global optical model po- 
tential of Koning and Delaroche 15]. 

For each event generated, FREYA first considers the 
possibility of pre-equilibrium neutron emission and, if 
it occurs, a neutron is emitted with an energy se- 
lected from the calculated pre-equilibrium spectrum (see 
Fig. [3]). Subsequently, the possibility of equilibrium neu- 
tron evaporation is considered, starting either from the 
originally agitated compound nucleus, ^^°Pu*, or the less 
excited nucleus, ^'^^Pu*, remaining after pre-equilibrium 
emission has occurred. Neutron evaporation is iterated 
until the excitation energy of a daughter nucleus is below 
the fission barrier (in which case the event is abandoned 
and a new event is generated) or the nucleus succeeds in 
fissioning. 



B. Mass and charge partition 

After the possible pre-fission processes, we are pre- 
sented with a fission- ready compound nucleus hav- 
ing an excitation energy Eq . The first task is to divide it 
into a heavy fragment Zh and a complementary light 
fragment "^^Z^. Since no quantitatively useful model is 
yet available for the calculation of the fission fragment 
mass yields, we have to invoke experimental data. We 
follow the procedure employed in the original version of 
FREYA [75. 

We thus assume that the mass yields Y{Af) of the 
fission fragments exhibit three distinct modes, each one 
being of Gaussian form [iGj . 



Y{Af) = SMf) + S2{Af) + S^Af) 



(4) 



The first two terms represent asymmetric fission modes 
associated with the spherical shell closure at iV = 82 
and the deformed shell closure at iV = 88, respectively, 
while the last term represents a broad symmetric mode, 
referred to as super- long [iTj . Although the symmet- 
ric mode is relatively insignificant at low excitations, its 
importance increases with the excitation energy and ul- 
timately dominates the mass distribution. 

The asymmetric modes have a two-Gaussian form. 



iAf-A+Dif/2af 



,(5) 



while the symmetric mode is given by a single Gaussian 
Nl 



Sl 



-iAf-Af/2al 



27r (Tz, 



(6) 



with A = ^Af). Since each event leads to two fragments, 
the yields are normalized so that — 2. Thus, 



2iVi + 2N2 + Nl 



(7) 



apart from a negligible correction because A/ is discrete 
and bounded from both below and above. 



Most measurements are of fission product yields [l8|, 
the yields after prompt neutron emission is complete. 
However, FREYA requires fission fragment yields, i.e. the 
probability of a given mass partition before neutron evap- 
oration has begun. While no such data are yet avail- 
able for Pu, there exist more detailed data for ^''^U(n, /) 
that give the fragment yields as functions of both mass 
and total kinetic energy, Y{Af,TKE) for E,, < 6 MeV 
p^ . Guided by the energy dependence of these data, 
together with other available data on the product yields 
from ^^^U(7i, /) |2Q,] and ^^^Pu(7i, /) 21] we derive an ap- 
proximate energy dependence of the fragment yields for 
239pu(n, /) up to En = 20 MeV. 

We now discuss how we obtained the values of the pa- 
rameters used in Eqs. ^ and ©. The displacements, 
Di, away from symmetric fission in Eq. ^ are anchored 
above the symmetry point by the spherical and deformed 
shell closures and because these occur at specific neutron 
numbers, we assume that the values of Di are energy 
independent. The fitted values of the displacements for 
235u(n, /) are Di = 23.05 and D2 = 16.54. The val- 
ues of Di should be smaller for ^^^Pu than for ^'^^U, 
Df — ZJf" 2, because the larger Pu mass is closer 
to the shell closure locations. We take Di = 20.05 and 
D2 — 14.54 for first-chance fission (Ao=240) and increase 
those values by i for each pre-fission neutron emitted. 

The widths of the asymmetric modes, tTj, are expanded 
to second order in energy: — + unEn + ai2E'^. 
We fix the energy dependence of Ci from the ^■^^U(n, /) 
fragment yields as a function of mass and total kinetic 
energy of Ref. [ll]. To adjust our results for ^^^U(n, /) 
to Pu, we assume that general energy dependence of the 
parameters is the same even though the values at En = 
are different. We find 

(71 = 5.6 + 0.0937 (£;„/MeV) + 0.034 (K/MeV)^ (8) 
(72 = 2.5 + 0.11060 (£;„/MeV) 4- 0.008 (i;„/MeV)^ (9) 

When the fissioning nucleus is the original system, ^^°Pu, 
then En is the value of the actual incident neutron en- 
ergy. But when the fissioning nucleus is lighter, i.e. when 
1^0 pre-fission neutrons have been emitted, then En is the 
equivalent incident neutron energy, i.e. the incident en- 
ergy that would generate the given excitation energy Eq 
when absorbed by the nucleus ^•^^"""Pu. The width of 
the super- long component, ctl, is assumed to be constant, 
independent of both the incident energy and the fission- 
ing isotope. We take — 12. 

The normalizations Ni change only slowly with inci- 
dent energy until symmetric fission becomes more prob- 
able, after which they decrease rapidly. We therefore 
model the energy dependence of Ni by a Fermi distribu- 
tion, 



N, ^N^{1 + exp[(S„ - Ei)/E2]y 



(10) 



We assume that the midpoint and the width are the same 
for both modes, Ei = 10.14 MeV and E2 = 1.15 MeV, 
so that the relative normalizations for the asymmetric 
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FIG. 4: (Color online) Calculated fragment yields as function 
of fragment mass for thermal (solid) and 14 MeV (dashed) 
neutrons. The 14 MeV result also includes contributions from 
multichance fission events. 
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FIG. 5: (Color online) The measured average TKE as a func- 
tion of the mass number of the heavy fragment [2^-[13 com- 
pared to FREYA calculations at thermal energies. The FREYA 
result is shown with the calculated variance around each A. 



modes have the same energy dependence. We have not 
assumed that E\ and are identical for U and Pu, be- 
cause the transition from asymmetric to more symmet- 
ric fission is not as smooth a function of energy in the 
few-MeV region for Pu as it is for U [13, [2l| . We take 
iVf = 0.757 and = 0.242. With TVi and N2 given by 
Eq. ([To]), Nl is determined from Eq. (0). 

The resulting fragment yields for two representative 
incident neutron energies are shown in Fig. 21 The deep 
dip at ^Aq visible for the thermal yields has substantially 
filled in by = 14 MeV. The fragment yield at 14 MeV 
is a composite distribution because there are substan- 
tial contributions from second- and third-chance fission 
for incident neutrons of such high energy (see Fig. [ij. 
The dashed curve in Fig. 2] includes these contributions 
weighted with the appropriate probabilities shown in 

Fig.m 

The overall broadening of the yields is due in part to 
the larger widths of the Si and S2 modes at higher ener- 
gies and in part to the increased contribution of the Sl 
component. We note that while ctl does not change, the 
larger enhances the importance of this component. 

Once the Gaussian fit has been performed, it is 
straightforward to make a statistical selection of a frag- 
ment mass number A f . The mass number of the partner 
fragment is then readily determined since A^ + Ah — 
Ao - 1^0- 

The fragment charge, Zf, is selected subsequently. For 
this we follow Ref. [23] and employ a Gaussian form, 

PA^Zf) oc e-(^/-^/(^/»'/2-l , (11) 

with the condition that \Zf — Zf{Af)\ < 5az- The cen- 
troid is determined by requiring that the fragments have, 
on average, the same charge-to-mass ratio as the fission- 
ing nucleus, Zf{Af) = AfZo/Ao. The dispersion is the 
measured value, az — 0.5 23]. The charge of the com- 
plementary fragment then follows using Zl + Zh = Zq. 



C. Fragment energies 

Once the partition of the total mass and charge among 
the two fragments has been selected, the Q value asso- 
ciated with that particular fission channel follows as the 
difference between the total mass of the fissioning nucleus 
and the ground-state masses of the two fragments, 

Qlh = Mi^^"-''°Fn*) - Ml - Mh . (12) 

FREYA takes the required nuclear ground-state masses 
from the compilation by Audi et al. |23l] . supp lemented 
by the calculated masses of Moller et al. [2J| when no 
data are available. The Qlh value for the selected fis- 
sion channel is then divided up between the total kinetic 
energy (TKE) and the total excitation energy (TXE) of 
the two fragments. The specific procedure employed is 
described below. 

Through energy conservation, the total fragment ki- 
netic energy TKE is intimately related to the resulting 
combined multiplicity of evaporated neutrons, i^l + i^H, 
which needs to be obtained very accurately. 

Figure [5] shows the measured average TKE value as 
a function of the mass number of the heavy fragment, 
Ah ■ Near symmetry, the plutonium fission fragments are 
mid-shell nuclei subject to strong deformations. Thus the 
scission configuration will contain significant deformation 
energy and a correspondingly low TKE. At Ah = 132, 
the heavy fragment is close to the doubly-magic closed 
shell having Zh = 50 and Nh = 82 and is therefore resis- 
tant to distortions away from sphericity. Consequently, 
the scission configuration is fairly compact, causing the 
TKE to exhibit a maximum even though the complemen- 
tary light fragment is far from a closed shell and hence 
significantly deformed. 

The TKE values shown in Fig. [5] were obtained in 
experiments using thermal neutrons (25l-l27|. Unfortu- 
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FIG. 6: (Color online) The average fragment kinetic energy 
as a function of fragment mass from Refs. [2^ . [27i] compared 
to FREYA calculations at thermal energies. The FREYA result 
is shown with the calculated variance around each A. 



nately, there are no such data for higher incident energy. 
We therefore assume the energy-dependent average TKE 
values take the form 



TKE{AH,En) = TKEdata(AH)+dTKE(£;„) . (13) 

The first term is extracted from the data shown in Fig. [SJ 
while the second term is a parameter adjusted to ensure 
reproduction of the measured energy- dependent average 
neutron multiplicity, V{En). In each particular event, the 
actual TKE value is then obtained by adding a thermal 
fluctuation to the above average, as explained later. 

Figure [5] includes the average TKE values calculated 
with FREYA at thermal energies, together with the as- 
sociated dispersions (these bars are not sampling errors 
but indicate the actual width of the TKE distribution for 
each Ah)- 

Figure [6] shows the single- fragment kinetic energy ob- 
tained with FREYA for incident thermal neutrons. Al- 
though FREYA is not tuned to match the single fragment- 
kinetic energies, it does reproduce these data quite well. 
Due to momentum conservation, the light fragment car- 
ries away significantly more kinetic energy than the heavy 
fragment. Furthermore, the kinetic energy of the frag- 
ment is nearly constant for Af < 106, but after the dip 
near symmetry there is an approximately linear decrease 
in the fragment kinetic energy. The figure also shows the 
calculated width in the fragment energy distribution, to- 
gether with a few typical experimental widths provided 
by Ref. [261]. 

Of particular interest is the dependence of the average 
neutron multiplicity on the fragment mass number Af, 
shown in Fig. [71 It is seen that the FREYA calculations 
provide a rather good representation of the 'sawtooth' be- 
havior oi D{Af), as shown in Fig. [71 even though FREYA is 
also not tuned to these data. Although the agreement is 
good, the observed behavior is not perfectly reproduced. 
When Af > 150, a region where the fragment yield is 
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FIG. 7; (Color online) The measured avera ge n eutron mul- 
tiplicity as a function of the fragment mass [26l - [2a ] together 
with thermal FREYA results for a = A/eo (circles) and from 
Eq. (|15|) . with both the average and the dispersion indicated 
for each Ai. 



decreasing sharply, the data and the calculations appear 
to diverge. We note that the uncertainties on the data 
in this region, where reported, are rather large. 

Once the average total fragment kinetic energy has 
been obtained, the average combined excitation energy 
in the two fragments follows by energy conservation. 



TXE 



}lh - TKE 



(14) 



The first relation indicates that the total excitation en- 
ergy is partitioned between the two fragments. As is 
common, we assume that the fragment level densities are 
of the form pi{E*) ^ exp{2\/ aiUi) , where Ui is the effec- 
tive statistical energy in the fragment. 

We have studied two forms of the level-density param- 
eter, Qi. One is the simple form ai = Ai/eo used in the 
Madland-Nix model; it has no energy dependence and 
may be appropriate for those (far-from-stability) frag- 
ments for which little information is available, even for 
the ground states. As an alternative, we have also used 
a parameterization based on the back-shifted Fermi gas 
(BSFG) model [H, 



Co 
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--rUi 



(15) 



where Ui = E* — Ai and 7 — 0.05, also used in Ref. [30|. 
The pairing energy of the fragment, , and its shell cor- 
rection, SWi, are tabulated in Ref. [23 based on the mass 
formula of Koura et al. fsijl . If SWi is negligible or if U 
is large then the renormalization of Ai is immaterial and 
the BSFG level-density parameter reverts to the simple 
form, Qi « Ai/eo. [Because the back shift causes Ui to 
become negative when E* is smaller than A^, we replace 
Ui{E*) by a quadratic spline for E* < 2Ai while retain- 
ing the expressions Ti = ^/uTJoi for the temperature and 
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cr|;. = 2UiTi for the variance of the energy distribution 
to ensure a physicaUy reasonable behavior.] In both sce- 
narios, we take cq as an adjustable model parameter. 

If the two fragments are in mutual thermal equilibrium, 
Tl — Th, the total excitation energy will, on average, be 
partitioned in proportion to the respective heat capaci- 
ties which in turn are proportional to the level-density 
parameters, i.e. ^ Qi. FREYA therefore first assigns 
tentative average excitations based on such an equipar- 
tition, 



E* 



aL{El)+aH{E'h) 



■TXE 



(16) 



where E* = {Ai/Ao)TXE. Subsequently, because the ob- 
served neutron multiplicities suggest that the light frag- 
ments tends to be disproportionately excited, the average 
values are adjusted in favor of the light fragment. 



El = xEl 



Eh = TKE - E^ 



(17) 



where x is an adjustable model parameter expected 
be larger than unity, as su gges ted by measurements of 
235U(n,f) m and 252cf(sf) |3l. 

After the mean excitation energies have been assigned, 
FREYA considers the effect of thermal fluctuations. The 
fragment temperature Ti is obtained from Ui = Ui{E*) = 
aiTf and the associated variance in the excitation E* is 
taken as a, = 2UiT,, where U{E*) = E* in the simple 
(unshifted) scenario. 

Therefore, for each of the two fragments, we sample a 
thermal energy fluctuation SE* from a normal distribu- 
tion of variance af and modify the fragment excitations 
accordingly, arriving at 



E* = E^+SE* , t = L,H. 



(18) 



Due to energy conservation, there is a compensating op- 
posite fluctuation in the total kinetic energy, so that 



TKE 



TKE - SE*r 



5E 



H 



(19) 



With both the excitations and the kinetic energies of 
the two fragments fully determined, it is an elementary 
matter to calculate the magnitude of their momenta with 
which they emerge after having been fully accelerated by 
their mutual Coulomb repulsion [3] . The fission direction 
is assumed to be isotropic (i.e. directionally random) in 
the frame of the fissioning nucleus and the resulting frag- 
ment velocities are finally Lorentz boosted into the rest 
frame of the original 240p^* gyg^gm. 



D. Fragment de-excitation 



[sot foi' ^^^Cf (sf) and 2^^U(n, /). After neutron emission 
is no longer energetically possible, FREYA simulates the 
sequential emission of photons by a similar method Q, 
see also Ref. 1341. 



1. Neutron evaporation 

Neutron emission is treated by iterating a simple neu- 
tron evaporation procedure for each of the two fragments 
separately. At each step in the evaporation chain, the 
excited mother nucleus "^^Zi has a total mass equal to 
its ground-state mass plus its excitation energy, M* — 
Mf^ + E* . The Q- value for neutron emission from the 
fragment is then Qn = M* — Mf — nin, where Alf is 
the ground-state mass of the daughter nucleus and to„ 
is the mass of the neutron. (For neutron emission we 
have Af — Ai — 1 and Zf = Zi) The Q-value is equal to 
the maximum possible excitation energy of the daughter 
nucleus, achieved if the final relative kinetic energy van- 
ishes. The temperature in the daughter fragment is then 
maximized at T™^^. Thus, once Qn is known, the ki- 
netic energy of the evaporated neutron may be sampled. 
FREYA assumes that the angular distribution is isotropic 
in the rest frame of the mother nucleus and uses a stan- 
dard spectral shape HH , 



UE) 



1 dN„ 



Ee-^'/^r 



(20) 



Nn dE 

which can be sampled efficiently Q- 

Although relativistic effects are very small for neutron 
evaporation, they are taken into account to ensure exact 
conservation of energy and momentum, which is conve- 
nient for code verification purposes. We therefore take 
the sampled energy E to represent the total kinetic en- 
ergy in the rest frame of the mother nucleus, i.e. it is 
the kinetic energy of the emitted neutron plus the re- 
coil energy of the excited residual daughter nucleus. The 
daughter excitation is then given by 



E*, 



Qr. 



E . 



(21) 

E^. The mag- 



and its total mass is thus — Mf" 
nitude of the momenta of the excited daughter and the 
emitted neutron can then be determined [7]. Sampling 
the direction of their relative motion isotropically, we 
thus obtain the two final momenta which are subse- 
quently boosted into the overall reference frame by the 
appropriate Lorentz transformation. 

This procedure is repeated until no further neutron 
emission is energetically possible, i.e. when E^^ < 5„, 
where Sn is the neutron separation energy in the prospec- 
tive daughter nucleus, Sn = M{^'' Zd)-M{^''-^ Zd)-mn. 



Usually both fully accelerated fission fragments are ex- 
cited sufficiently to permit the emission of one or more 
neutrons. We simulate the evaporation chain in a man- 
ner conceptually similar to the method of Lemaire et al. 



2. Photon radiation 

After the neutron evaporation has ceased, the ex- 
cited product nucleus may de-excite by sequential photon 
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emission. FREYA treats this process in a manner analo- 
gous to neutron evaporation, i.e. as the statistical emis- 
sion of massless particles. While this simple treatment is 
expected to be fairly reasonable at the early stage where 
the level density can be regarded as continuous, it would 
obviously be inadequate for late stages that involve tran- 
sitions between specific levels. 

There are two important technical differences relative 
to the treatment of neutron emission. There is no sep- 
aration energy for photons and, since they are massless, 
there is no natural end to the photon emission chain. 
We therefore introduce an infrared cutoff of 200 keV. 
Whereas the neutrons may be treated by nonrelativistic 
kinematics, the photons are ultrarelativistic. As a conse- 
quence, their phase space has an extra energy factor. 



ME) 



1 dN-f 
'N^~dE 



(22) 



The photons are assumed to be emitted isotropically 
and their energy can be sampled very quickly from the 
above photon energy spectrum Q. The procedure is 
repeated until the available energy is below the spec- 
ified cutoff, yielding a number of kinematically fully- 
characterized photons for each of the product nuclei. 



III. RESULTS 

We now proceed to discuss our analysis of the prompt 
fission neutron spectrum (PFNS). We first describe the 
computational approach and then explain how the model 
parameters are determined. The resulting prompt neu- 
tron spectral evaluations are then discussed in detail. Fi- 
nally, we present some additional observables of particu- 
lar relevance. 



A. Computational approach 

Here we briefly describe the statistical method used for 
determining model parameters and reaction observables. 
Our analysis uses the Monte Carlo approach to Bayesian 
inference outlined in many books on general inverse prob- 
lem theory, e.g. Ref j36i] . 

We have introduced several model parameters: eo, x, 
and dTKE, which in principle may be adjusted for each 
incident neutron energy i?„. However, because indepen- 
dent fits to the experimental data tend to yield values of 
eo and x that are nearly independent of we shall 

assume that these two parameters are energy indepen- 
dent. This simplification will facilitate the optimization 
procedure. Thus a given model realization is character- 
ized by the two values eo and x together with the function 
dTKF,{En) which, for practical purposes, will be defined 
by its values at certain selected energies, {dTKE^}. For 
formal convenience, we denote the set of model parame- 
ter values as m = {mk}. 



When FREYA is used with any particular value set m, 
it yields a sample of fission events from which we can ex- 
tract observables, d{m), that can be directly compared 
to the corresponding experimental values, tZexp- For ex- 
ample, we may extract the energy-dependent mean neu- 
tron multiplicity, F(£'„), and compare it with the values 
given in the ENDF/B-VH.O evaluation 0. 

We assume that the experiment provides not just the 
values but also the entire associated covariance matrix 
Sexp- (The diagonal elements of Scxp are the variances 
on the individual observables.) The degree to which the 
particular model realization defined by the parameter 
values m describes the measured data d^^p is then ex- 
pressed by 

(23) 



P(<ioxp|m) - exp {~^x^{m)) 



where x^(^) is the generalized least-squares deviation 
between the model m and experiment, 

= (dcxp - d{m)) ■ (Sexp)"' • (dcxp - d{m)f . (24) 
Employing merely the diagonal part of S 



oxpj 



the 



certainties alone, ensures that well-measured observables 
carry more weight than poorly measured ones. This ap- 
proach was used in the previous PFNS evaluation [&\, 
which was restricted to lower energy (£"„ < 5.5 MeV). 
Here we now employ the full covariance matrix, thereby 
ensuring that correlations between measured observables 
are also taken into account. As we shall see, these corre- 
lations do impact the results. 

Using the above framework, we may now compute the 
weighted averages of arbitrary observables O — {Oi}. 
We assume that the physically reasonable values of the 
model parameters m are uniformly distributed within a 
hypercube in parameter space. This defines the a-priori 
model probability distribution P(rn). The best estimate 
of the observable Oi is then given by 



^ ^ = J dmP(m)P(dexp|m)C',(m) 



(25) 



The best estimate for the covariance between two such 
observables can be obtained similarly, 

EE ^ 0,0j > <0,>--<Ojy 

^ ^0^- ^ O, Oj- -< Oj . (26) 

In particular, we may compute the best estimate of V, 



^ v y = J dm P{rn) P{dcxp\rn) iy{rn) 



(27) 



and the prompt neutron spectrum, as well as the covari- 
ances between those quantities. 

In practice, we average over parameter space employ- 
ing a Monte Carlo approach, thereby reducing the inte- 
gral over all possible parameter values m to a sum over 
N sampled model realizations, {m'^")}. 



1 ^ 

^O,^^ -^P(m("))P(dexp|m("))0,(m(")) 



(28) 
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The joint probability w„ = P(m'^"^)P(<icxp|'Ti(")) may 
be viewed as the hkehhood that the particular model 
realization m'") is "correct." Since it depends exponen- 
tially on the likelihood tends to be strongly peaked 
around the favored set. It is important that the parame- 
ter sample be sufficiently dense in the peak region to en- 
sure that many sets have non-negligible weights. We use 
Latin Hypercube sampling (LHS) [33, IH], which samples 
a function of K variables with the range of each variable 
divided into M equally-spaced intervals. Each combina- 
tion of M and K is sampled at most once, with a maxi- 
mum number of combinations being {Ml)^~^. The LHS 
method generates samples that better reflect the distri- 
bution than a purely random sampling would. Conse- 
quently, relative to a simple Monte Carlo sampling, the 
employed sampling method requires fewer realizations to 
determine the optimal parameter set. We used 5000 re- 
alizations of the parameter space to obtain the optimal 
parameter values. 

Even though both V and the neutron spectrum are im- 
portant observables, the fact that the uncertainties on the 
evaluated V are so relatively small drive the results. In- 
deed, we use the evaluated V to constrain our new evalu- 
ation of the PENS, as in Ref. f6|. Thus, in our treatment, 
the spectra is an outcome rather than a comparative ob- 
servable. We use the evaluated V in the ENDF/B-VII.O 
database Q with the covariance resulting from the least- 
squares fit to the available ^^^Pu(n, /) data described 
in Ref. [s^. The energy- dependent neutron multiplicity, 
77(i5„) is represented as a locally linear fit to the experi- 
mental data. Since the nodes in this fit do not align with 
fitted data the fit introduces energy correlations that are 
encoded into the covariance matrix. 

Eor each model realization, FREYA is used to generate a 
large sample of fission events (typically one million events 
for each parameter set) for each of the selected incident 
neutron energies and the resulting average multiplicity 
77(£'„) is extracted from the generated event sample. 



B. Fit results 

Given the tendency of cq to be larger than 8 MeV, 
we let Co vary between 8 and 12 MeV. Also, recalling our 
previous results ^] and the experimental indications that 
the light fragment is hotter than the heavy fragment, we 
have assumed 1 < a: < 1.4. The resulting optimal val- 
ues of these parameters are eo = 10.0724 ± 0.5571 MeV 
and X = 1.2339 ± 0.0496 for the BSFG level density and 
eo = 9.6068 ± 0.8930 MeV and x = 1.0164 ± 0.1661 for 
a = A/ Co- These values of eo are consistent with the 
calculation of a in Eq. P3|) which does not include col- 
lective effects. If collective behavior is included, then we 
expect eo ^ 13 MeV based on the RIPL-3 systematics 
[40l]. Previously, we obtained eo 8 MeV and x ~ 1.1 
with a slightly different TKE prescription and usiiig only 
the diagonal elements of the V covariance matrix Q . 

We have fit dTKE at eight values of incident neutron 
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FIG. 8; (Color online) The fitted values of dTKE as a func- 
tion of incident neutron energy for a — A/eo (circles) and 
the BSFG (squares). The locations of the node points are 
indicated by diamonds at dTKE — 1. 

energy, = 10"", 0.25, 2, 4, 7.5, 10.5, 14 and 20 MeV, 
to keep the parameter space manageable. These points 
are chosen to reflect the physics of the fission process: the 
region between 0.25 and 2 MeV is where v{En) changes 
slope while the second and third-chance fission thresholds 
occur in the regions 4 < i?„ < 7.5 MeV and 10.5 < 
En < 14 MeV, respectively. The full 20-point grid of the 
FREYA evaluation is then covered by means of a linear 
interpolation between these node points. The resulting 
values of dTKE for the two level density prescriptions 
are shown in Fig. [S] The locations of the node points 
are indicated by diamonds at dTKE = 1. The error bars 
on dTKE at the node points are the standard deviations 
obtained from the averaging over the range of parameter 
values while the error bars on dTKE between two node 
points are the interpolated dispersions between those two 
points. The standard deviations of dTKE are larger for 
the case where a = A/e^^ possibly because the value of 
X is close to the edge of the fit range. 

Because dTKE represents the shift in the total frag- 
ment kinetic energy from the value obtained for incident 
thermal neutron energies, dTKE should depend on the 
incident neutron energy, as suggested in Eq. (1131) . The 
value of dTKE is positive, indicating that using the ther- 
mal average value of TKE leads to too many neutrons. 
The positive dTKE is then required to reduce the excita- 
tion energy sufficiently to give a good fit to V. For exam- 
ple, reducing dTKE at E^ = 0.5 MeV from 1.4 MeV to 
zero while retaining the same values of eo and x, reduces 
TKE by about 1% and increases 77 by ^ 6%. The change 
in dTKE with energy is not particularly systematic in 
either scenario. 

Above 14 MeV, the ENDF/B-VII.O v evaluation is not 
based on data but on a linear extrapolation of measure- 
ments taken for higher incident energies. Thus, v is not 
well constrained near the high end of the energy range. 

A direct comparison of our fitted values of V with those 
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FIG. 9: (Color online) The difference between the ENDF- 
B/VII.O evaluation and our fits to V using FREYA with a — 
A/cq (circles) and the BSFG result (squares). The locations 
of the node points are indicated by diamonds at (!^)endf — 

(i^)FREYA = 0. 



in the ENDF-B/VII.O evaluation is less revealing than 
than fit residuals, the difference i^endf ~ i'freya- The 
residual values are shown in Fig. [HI The standard de- 
viation on each point reflects only the uncertainty on 
the ENDF-B/VII.O evaluation and not the uncertainty 
on the fitted V. The uncertainties on V arising from the 
fitting procedure are typically smaller than those on the 
evaluated values. The largest residuals occur for En > 2 
MeV. The large uncertainties on the residuals at 16 and 
20 MeV arise from the lack of experimental data at these 
values of E^- 

We note that the ENDF-B/VII.O V evaluation lies more 
than one standard deviation above the evaluated v data 
extracted in the ENDF-B/VII.O covariance analysis in 
the region 0.1 < £;„ < 1 MeV Our results with both 
calculations of the level-density parameter agree rather 
well with the evaluated data in this region. Here, where 
j7 is smallest, the relative difference is less than 0.5%. At 
higher energies the relative uncertainty increases because, 
while V increases with incident energy, the residual also 
increases. The residual difference is largest for the BSFG 
at 6 < £^n < 14 MeV. We note that, assuming each 
energy point is independent of all others, e.g. the errors 
are uncorrelated and the covariance matrix is diagonal, 
then the fit residuals can be very small, as in our previous 
work 0. 

Examples of the resulting prompt fission neutron spec- 
trum for the BSFG level density are shown in Fig.[TUl We 
present dv/dE for four representative energies: En =0.5 
MeV (thermal neutrons), 4 MeV (below the second- 
chance fission threshold), 9 MeV (below the threshold 
for third-chance fission), and 14 MeV (relevant for certain 
experimental tests). We note that the integral of dv/dE 
over outgoing neutron energies gives the average neutron 
multiplicity, 17, obtained from the fitting procedure. In 
addition to the increase in the peak of the spectrum, we 
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FIG. 10: (Color online) The prompt neutron spectra at se- 
lected energies resulting from our BSFG fits. 



note that the average outgoing neutron energy appears 
to increase with incident energy. 

Figure [Til compares results for the BSFG level density 
prescription and a — A/cq at En = 0.5 and 14 MeV. 
While the peaks appear at the same point at 0.5 MeV, 
the spectra with a = A/cq is broader around the peak 
than the BSFG case because the back shift has the effect 
of narrowing the spectrum. This effect is particularly 
apparent at 14 MeV where the BSFG peak is significantly 
higher than that for a = A/cq even though the residual 
differences between the evaluation and the fitted 77 values 
in Fig. iniare negligible on this scale. 

A close inspection of the spectra obtained for 9 and 
14 MeV will reveal abrupt drops in value at the ener- 
gies corresponding to the threshold for emission of a sec- 
ond pre-fission neutron, namely at E2 = En — S'n(^^^Pu), 
3.4 and 8.4 MeV, respectively. When the energy of the 
first pre-fission neutron is below E2, the daughter nu- 
cleus is sufficiently excited to make secondary emission 
possible. (These threshold discontinuities are also visi- 
ble in the spectral differences shown in Figs. \TI\ and [T51) 
This effect, noted already by Kawano et al. |10||, grows 
larger at higher incident energies where multichance fis- 
sion is more probable. Furthermore, when the com- 
bined energy of the first two pre-fission neutrons is below 
E3 = En- S'„(239pu) - S^i^^^Pu) the emission of a third 
pre-fission neutron is energetically possible. The thresh- 
old £^3 appears as a change in the FREYA slope relative 
to the continuous neutron spectrum in ENDF-B/VII.O, 
as is visible near 1.35 MeV in the En=14: MeV curve in 
Figs. [T^ and [T^ In principle, these onset effects can be 
measured experimentally which could thus provide novel 
quantitative information on the degree of multichance 
fission. 

Our final evaluation is based on our fits to V and its 
energy-energy covariance, as described above. The re- 
sulting spectral evaluations for both level density scenar- 
ios are incorporated in different ENDF-type files with 
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FIG. 11: (Color online) Comparison of the spectra resulting 
from the BSFG and a = A/eo level density parameterizations 
at En = 0.5 and 14 MeV. 



their spectral shapes alone. 

To produce the spectral evaluation, requiring fine en- 
ergy spacings over the range 10~^ < £^ < 20 MeV, from 
our FREYA results, we fit the FREYA PENS in the regions 
where either the employed bin widths are not sufficiently 
small (in the region below 0.1 MeV) or the statistics are 
limited (above 6-8 MeV). In between, we interpolate the 
FREYA spectra. We also interpolate the spectra around 
the multichance fission thresholds where smooth fitting 
would not be appropriate. 

Figures [T^ and [T3] give the difference between the 
present evaluations and the ENDF/B-VII.O evaluation. 
Both are normalized to unity at each value of En- For 
incident energies below the threshold for multichance fis- 
sion, the difference between the two evaluations is around 
10% for both methods of calculating the level density. 
In the case where a — A/cq, our spectra are systemat- 
ically higher than the ENDF/B-VII.O result. With the 
BSFG parameterization of the level density, the FREYA 
spectra are systematically lower than the ENDF-B/VII.O 
evaluation. As the incident energy increases up to ~ 14 
MeV, in both cases the FREYA spectra are generally lower 
than ENDF/B-VII.O for outgoing energies between 0.01 
and 3 MeV because pre-fission neutron emission lowers 
the effective temperature of the fragments when fission 
occurs. The FREYA evaluation has a higher-energy tail 
than ENDF/B-Vll.O for all E,,, giving the FREYA resuhs 
a higher average energy, particularly for the BSFG. 

As noted in the discussion of Figs. [TU] and [Til con- 
tributions from pre-fission neutron emission change the 
calculated spectral slope a.t E = En — Sn- Although 



the Madland-Nix (Los Alamos model) evaluation [11[ in- 
cludes an average result for multichance fission, the spec- 
tral shape for pre-fission emission is assumed to be the 
same as that of prompt neutron emission (evaporation) 
post fission. Thus the ENDF/B-VII.O evaluations are 
always smooth over the entire outgoing energy regime, 
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FIG. 12: (Color online) The percent difference between the 
ENDF-B/VII.O evaluation of the PFNS and our BSFG result 
at several representative incident neutron energies. Results 
are shown for 0.5 (solid black), 4 (short dashed red), 6 (dot- 
dashed blue), 10 (dot-dot-dashed green), 14 (dot-dash-dashed 
violet) and 20 (long dashed magenta) MeV. 
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FIG. 13: (Color online) The percent difference between the 
ENDF-B/VII.O evaluation of the PFNS and our result with 
a = j4/eo at several representative incident neutron energies. 
Results are shown for 0.5 (solid black), 4 (short dashed red), 
6 (dot-dashed blue), 10 (dot-dot-dashed green), 14 (dot-dash- 
dashed violet) and 20 (long dashed magenta) MeV. 



regardless of incident energy, while the FREYA evalua- 
tions refiect the changes due to pre-fission emission. The 
slope changes at the multichance fission thresholds in the 
FREYA spectra are evident for the i5„ = 6, 10 and 14 MeV 
difference curves at 0.35, 4.35 and 8.35 MeV, respectively. 
We note that the threshold at 0.35 MeV is not well de- 
fined for a = A/tQ while it is exaggerated for the BSFG. 

FigureHHshows the ratio of our 0.5 MeV results for the 
a = A/eo and BSFG scenarios to a Maxwell distribution 
with T = 1.42 MeV. The ENDF-B/VII.O ratio is also 
shown, along with data from Refs. pl| - |49| . As expected 
from the comparison in Figs. [T^] and [T31 the FREYA 
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FIG. 14: (Color online) The ratios of ENDF-B/VII.O (solid), 
FREYA BSFG (long dashed), and FREYA a = A/eo 
(dot-dashed) relative to a Maxwell distribution with T — 
1.42 MeV. The data from Refs. [ill-ligl] are also shown as 
ratios to the same Maxwellian. 

ratios bracket the low energy ENDF-B/VII.O ratio. At 
E > 0.3 MeV, the dip relative to ENDF-B/VII.O around 
1 MeV followed by a bump at ^ 3 MeV and subsequent 
decrease in the percent difference in Fig. [T^ for the BSFG 
is seen reflected in the ratio relative to the Maxwellian. 
At higher outgoing energies, the BSFG result is higher 
than the Maxwellian. On the other hand, the a = A/eq 
result is similar to that of ENDF-B/VII.O but with a 
slightly lower average energy. Since there is a great deal 
of scatter in the data over the entire energy range, any 
conclusion about the quality of the fits with respect to 
the spectral data is difficult. 

C. Covariances 

We can calculate covariances and correlation coeffi- 
cients between the optimal model parameter values as 
well as between the various output quantities using Eq. 
(PS)) . The covariance between two parameters ruk and 
mk' is 

Sfcfc' = -< {nik < mk >-){mk' < ruk' >-) > ■ (29) 

The diagonal elements, Yiku are the variances S^^^ , rep- 
resenting the squares of the uncertainty on the optimal 
value of the individual model parameter mfc, while the 
off-diagonal elements give the covariances between two 
different model parameters. It is often more instruc- 
tive to employ the associated correlation coefficients^ 
Ckk' = Sfcfc'/PnifcSmfe']i which is plus (minus) one for 
fully (anti)correlated variables and vanishes for entirely 
independent variables. The correlation coefficient be- 
tween the energy independent parameters eo and x is 
Ceo, a; ~ 0-985, a near-perfect correlation, for a — A/cq, 
while C'eo.x ~ —0.7, a strong anticorrelation, for the 
BSFG level density. 



FIG. 15: (Color online) The correlation coefficients between 
dTKE(i5„) and either eo (circles) or x (triangles) for a = A/eo 
and eo (squares) or x (diamonds) for the BSFG level density. 



The correlation coefficients CdTKE,x and CdTKE,eo are 
shown in Fig. [15] When a — A/eo, the parameters are 
all strongly correlated although CdTKE,ea is slightly larger 
than CdTKE,x- On the other hand, the BSFG inputs are 
either strongly correlated or anticorrelated with the sign 
of the correlation changing near the second and third- 
chance fission thresholds. The anticorrelation of CdTKE.x 
and CdTKE.eo with each other is consistent with the an- 
ticorrelation in Ceo.x mentioned previously and may be 
related to the larger value of x in this case. 

We may also compute the covariance between the 
spectral strength at different outgoing energies E using 
Eq. (pS)) . The resulting correlation coefficients Cei,E2 
are shown in Fig. [TBI and [T7I for En = 0.5 MeV. We see 
that Cex,E2 ~ 1 (light areas) when the two specified en- 
ergies lie on the same side of the crossover region, while 
Cei,E2 ~ — 1 (dark areas) when they lie on opposite sides. 
The crossover region is from 2.5 to 4 MeV indicates that 
the spectrum tends to pivot around E fa 3.5 MeV when 
the parameter space is explored, slightly higher E than 
in Ref. [|. 

The BSFG correlation in Fig.[Tn]is rather noisy. In con- 
trast, the crossover for a — A/eo in Fig. [T7|is very sharp 
with threshold-like boundaries between regions of almost 
perfect correlation and anticorrelation. In our previous 
paper @ , we used Eq. ([T5|) with U = E* , thus the results 
shown there exhibited slightly weaker correlations than 
shown in Fig. [TTI 

Figures \TE\ and [T9] show cuts at constant total neutron 
energy, Ek + Ek' ■ Similar results are found for all other 
incident energies considered. The trends are the same 
for both scenarios, large positive correlations at low AE, 
becoming negative at larger AE. This is because when 
model parameters are varied, the spectral shapes pivot 
about a single energy, Spivot ^3.5 MeV in these calcu- 
lations. Thus when both Ei and E2 are less than E'pivot, 
the differential changes are in phase and Cei.E2 ^ 1- If 
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Neutron energy El 

FIG. 16; (Color online) Contour plot of the correlation coeffi- 
cient (see Ea.l26|) between the spectral strengths at two differ- 
ent energies, Cei,E2, as obtained for En — 0.5 MeV with the 
BSFG result. The correlation changes from values near -1-1 
in the light regions (lower-left and central regions) to values 
near -1 in the dark regions (near the two axes). Each of the 
three straight lines connects points at which the two neutrons 
have the same combined energy, Ei + E2 = 5, 10, 15 MeV. 



e.g. El < iJpivot and E2 > -Epivot, differential changes in 
the spectra give an anticorrelation. 

The fluctuations in Cei,E2 are very large for the BSFG. 
There is a shift from negative back to positive and back 
to negative near the same value of AE where Cei,E2 
changes from -1-1 to —1 in Fig.[T21 This coincides with the 
behavior of the correlation matrix in Fig. [16] around the 
pivot point. The change in Cei.E2 is similar for Et — 5 
and 10 MeV, likely because the lines of constant Et sit 
at approximately equidistant locations on opposite sides 
of the pivot point. In our previous paper, these results 
were further apart because iJpivot ~ 2 MeV Q. 

There are several possible reasons why the BSFG level 
density causes fluctuations in Cei,E2- If a = A/cq, a 
rises smoothly with A and is energy independent. Thus, 
the temperature is also independent of incident energy. 
Strong correlations are also observed in other calculations 
based on average fission models such as Madland-Nix 
[5^ . Introducing the BSFG parameterization, Eq. (ITSI) . 
at fixed U , ignoring the pairing energy, as in our previ- 
ous paper, introduces fluctuations in a which soften the 
linear rise of a with A and reduce the sharpness of the 
correlations. Including the back shift due to the pairing 
energy further reduces the correlations, narrowing the 
peak in the spectra. Thus the back shift interferes with 
the correlations, introducing the noise shown in Fig. 1161 




Neutron energy E1 

FIG. 17: (Color online) Contour plot of the correlation coef- 
ficient (see Eq. I26() between the spectral strengths at two dif- 
ferent energies, Cex,E2, as obtained for En ~ 0.5 MeV with 
a = A/eo. The correlation changes rapidly from values near 
-1-1 in the light regions (lower-left and central regions) to val- 
ues near -1 in the dark regions (near the two ax;es). Each of the 
three straight lines connects points at which the two neutrons 
have the same combined energy, Ei + E2 = 5, 10, 15 MeV. 




Energy difference AE = |E^-E2l (MeV) 



FIG. 18: (Color online) The spectral correlation coefficients, 
Cei,E2, along the three lines of constant combined energy 
indicated in Fig. [16] for the BSFG result at En = 0.5 MeV. 



D. Benchmark tests 

There are several standard validation calculations that 
can be used to test our PFNS evaluations. They are criti- 
cal assemblies which test conditions under which a flssion 
chain reaction remains stationary, i.e. exactly critical; 
activation ratios which can be used to test the model- 
ing of the flux in a critical assembly; and pulsed sphere 
measurements which test the spectra for incident neutron 
energies of ~ 14 MeV. 

To perform these, we replace the ^^^^Pu PFNS in 
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FIG. 19: (Color online) The spectral correlation coefficients, 
Cei,E2, along the three lines of constant combined energy 
indicated in Fig. [l7] for the BSFG result at En = 0.5 MeV. 



the ENDL2011.0 database, identical to that in ENDF- 
B/VII.O, with our evaluations obtained with the BSFG 
level density parameterization and a = A/bq. In this sec- 
tion, we describe these tests and the FREYA results. As 
we will see, the two different treatments of a do about 
equally well on these benchmark tests with no conclusive 
preference. 



E. Validation against critical assemblies 

Critical assemblies, which are designed to determine 
the conditions under which a fission chain reaction is 
stationary, provide an important quality check on the 
spectral evaluations. The key measure of a critical as- 
sembly is the neutron multiplication factor kcs- When 
this quantity is unity, the assembly is exactly critical, 
i.e. the net number of neutrons does not change so that 
for every neutron generated, another is either absorbed 
or leaks out of the system. For a given assembly, the de- 
gree of criticality depends on the multiplicity of prompt 
neutrons, their spectral shape, and the energy-dependent 
cross section for neutron-induced fission. 

Plutonium criticality is especially sensitive to the 
prompt neutron spectrum because the ^'^^Pu(n, /) cross 
section rises sharply between En — 1.5 and 2 MeV, near 
the peak of the fission spectrum. As a result, increasing 
the relative number of low-energy neutrons tends to de- 
crease criticality, lowering k^g , while increasing the num- 
ber of higher energy neutrons increases criticality. 

Figure [201 shows calculations of kcs for four different 
plutonium assemblies from the criticality safety bench- 
mark handbook [5l|. Apart from the spectra, all data 
used in these calculations were taken from ENDF/B-VII 
Q. We show the kes for the two level density calcula- 
tions with our evaluated spectra and V from the ENDF- 
B/VII.O evaluation. The BSFG level density is super- 
critical (fcoff > 1) while the a — A/cq result is subcritical 
(fccff < !)• 

The energy-independent result of Ref. ^] led to values 
of fccff that were « 1.5 standard deviations away from the 
measured value for the Jezebel assembly which is sensi- 
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FIG. 20: (Color online) Calculated values of k^s for several 
^^^Pu critical assemblies obtained using our fits, labeled FREYA 
(solid squares and circles), in the Mercury code. The re- 
sults are compared to calculations using the ENDF-B/VII.O 
(blue diamonds) and proposed ENDF/B-VII. 1 (red squares) 
databases. 



live to fission induced by fast {En ~ 1 MeV) neutrons. 
By contrast, our results for both methods of calculating 
the level density, are within one standard deviation of the 
measured value and thus represent some improvement. 

We note that the other inputs to the Jezebel assembly 
test were highly tuned to match the k^s- The fact that 
replacing only the PENS without a full reevaluation of 
all the inputs to the criticality tests leads to a result 
that is inconsistent with fccff = 1 should therefore not 
be surprising. For example, if the average energy of the 
ENDF-B/VII.O PENS is high, the V evaluation would be 
forced higher to counter the effect on fceff. In addition, 
the inelastic (n, n') cross section is not well known and 
could require compensating changes that affect fcefi [s^ . 

To make further improvements in the evaluations with 
respect to assemblies, it would be useful to have an inline 
version of FREYA for use in the simulations. 
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F. Validation against activation ratios 

In the 1970's and 1980's, LANL performed a series of 
experiments using the spectra from the critical assem- 
blies Jezebel (mainly ^^^Pu), Godiva, BigTen and Flat- 
top25 (all primarily enriched uranium) to activate foils 
of various materials [11] . The isotopic content of the foils 
can be radiochemically assessed both before and after ir- 
radiation. Thus these measurements of the numbers of 
atoms produced per fission neutron are integral test of 
specific reaction cross sections in the foil material. Con- 
versely, a well-characterized material can also be used to 
test the critical assembly flux modeling. We have simu- 
lated several foils (239pu, ^^^\J, ^SSy, 238u^ 237]s^p^ Sly^ 

55Mn, 63Cu, 93Nb, Ag, i^isb, i39La, ^^^ir and ^^^Au) 
which are tests of the {n, /) and (n, 7) reactions in the 
Jezebel assembly. In all cases, our simulated values of 
the activation rates in these foils are either consistent 
with measured values or previous modeling using a mod- 
ified version of the ENDF/B-VII.O nuclear data library 
(5^ . As the fission cross sections for ^^^U, ^^^U, ^^^U, 
and ^^^Np are accurately known and span all incident 
neutron energies, these tests merely confirm our earlier 
modeling of plutonium critical assemblies. The neutron 
capture cross sections of the other foils are important for 
incident neutron energies less then 1 MeV but they are 
not known nearly as well as the fission cross sections. 
Thus our Calculated/Experiment ratios scatter around 
unity in these cases. 

Because both the Jezebel and Godiva critical assem- 
blies test the fast-neutron spectrum, with a significant 
portion of their neutron fluxes above 5 MeV, either might 
be used to test the high-energy portion of the ^^^Pu 
PFNS. Unfortunately, none of the tests that have been 
performed to date are useful for testing our FREYA PFNS 
evaluation. While many of the studied materials have 
high (> 10 MeV) (n, 2n) thresholds, the only (n, 2n) 
threshold material tested in Jezebel is ^^^Tm. Unfortu- 
nately, the "'^^^Tm (n, 2n) cross section is poorly known. 
There were also experiments with plutonium foils placed 
in uranium assemblies, but these tests are not useful for 
testing the ^^^Pu PFNS because the foils are too thin for 
secondary scatterings to play a significant role. It would 
particularly interesting to carry out a new set of (n, 2n) 
foil measurements using well-characterized materials in a 
primarily plutonium critical assembly to specifically test 
the high-energy portion of the spectrum. 



G. Validation with LLNL pulsed spheres 

The ENDL2011.0 database [H], including our FREYA 
evaluation, was tested against LLNL pulsed-sphere data, 
a set of fusion-shielding benchmarks [s^. The pulsed- 
sphere program, which ran from the 70's to the early 
90's, measured neutron time-of-flight (TOF) and gamma 
spectra resulting from emission of a 14 MeV neutron 
pulse produced by d-|-t reactions occurring inside spheres 
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FIG. 21: (Color online) The measured LLNL pulsed-sphere 
test data [531 (points) compared to calculations using either 
ENDF/B-VII.O or ENDL2011.0 with the present FREYA evaL 
nation included. 



composed of a variety of materials [57|. Models of 
the LLNL pulsed-sphere experiments using the Mercury 
Monte Carlo were dev elop ed for the materials reported 
in Goldberg et al. (ssl. l59l| . 

Figure compares results of our evaluation with 
the experimental data (STj and calculations based on 
ENDF/B-VII.O. The only difference between the Pu 
evaluations in these two calculations is the PFNS and 
associated F, all other quantities remain the same. Rel- 
ative to the ENDF/B-VII.O calculation, in both cases 
the FREYA spectra yields significantly better agreement 
with the data in the region around the minimum of the 
time-of-flight curve at 210 ns. However, it lies some- 
what higher than the ENDF/B-VII.O curve and the data 
during the subsequent rise (240 — 300 ns), particularly 
for the BSFG level densities. Thereafter, the two calcu- 
lations are practically identical and agree well with the 
data until w 480 ns. 

Such pulsed-sphere experiments test the PFNS at in- 
cident energies higher than those probed in critical as- 
sembly tests. The measured outgoing neutrons have a 
characteristic time-of-flight curve, see Fig.[5TJ The sharp 
peak at early times is due to 14 MeV neutrons going 
straight through the material without signiflcant inter- 
action. The dip at around 200 ns and the rise imme- 
diately afterward is caused by secondary scattering in 
the material as the neutrons travel out from the center. 
The location and depth of the dip is due to inelastic di- 
rect reactions, the high-energy tail of the prompt fission 
neutron spectrum, and pre-equilibrium neutron emission. 
The last two items are directly addressed by the present 
evaluation. A large part of the secondary interactions are 
due to neutrons that have interacted in the material and 
are thus less energetic than those from the initial 14 MeV 
pulse. At late times, where the results from both evalu- 
ations deviate from the data, the time-of-flight spectrum 
is dominated by scattering in surrounding material such 
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FIG. 22: (Color online) The total average excitation energy 
available for neutron emission (top), the average neutron ki- 
netic energy (middle) and the average neutron multiplicity 
(bottom) as a function of fragment mass number Af for ther- 
mal neutrons (black circles), En = 4 (red squares), 9 (green 
diamonds) and 14 (blue triangles) MeV. 



as detector components and concrete shielding. 



H. Additional observables 

The mass-averaged fragment kinetic energies obtained 
with FREYA are almost independent of the incident neu- 
tron energy E^. This feature is consistent with mea- 
surements made with ^ssy and ^38^ targets over similar 
ranges of incident neutron energy, 0.5<£'„<6 MeV [l9| 
and 1.2<i?„<5.8 MeV [6lj, respectively. In both cases, 
the average TKE changes less than 1 MeV over the entire 
energy range. 

However, Ref. [6l[ also showed that, while the mass- 
averaged TKE is consistent with near energy indepen- 
dence, higher-energy incident neutrons typically give less 
TKE to masses close to symmetric fission and somewhat 
more TKE for Ah > 140. Such detailed information is 
not available for neutrons on ^■^^Pu. We have therefore 
chosen to use a constant value of dTKE at each energy. 

Neutron observables are perhaps more useful for model 
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FIG. 23: (Color online) The total average excitation energy 
available for neutron emission (top), the average neutron ki- 
netic energy (middle) and the average neutron multiplicity 
(bottom) as a function of fragment mass number Af for in- 
cident neutrons at En = 0.5 MeV with the level densities 
calculated with a = A/eo (circles) and the BSFG approach 
(squares). The variances on the calculated results are shown 
for = 4 MeV only. 



validation. The near independence of TKE{Ah) on in- 
cident energy implies that the additional energy brought 
into the system by a more energetic neutron will be pri- 
marily converted into internal excitation energy. This is 
illustrated in the top panel of Fig. [22] which shows the 
fragment excitation E*{Af) for the BSFG level density 
with the representative incident energies considered in 
Fig. Hn] (thermal, 4, 9 and 14 MeV): E*{Af) increases 
linearly with incident energy. We note that the form of 
TKE{Ah) (see Fig. [5]) leads to the familiar sawtooth form 
ofE*{Af). 

The average kinetic energy of the evaporated neutrons 
is given hy E = 2T for a single emission, where T is 
the maximum temperature in the daughter nucleus, so 
T"^ (X E* — Sn for the first emission. Consequently, E 
should vary relatively little with Af, as is indeed borne 
out by the results for E{Af) shown in the middle panel 
of Fig. [221 The average outgoing neutron kinetic energy 
increases slowly with the incident neutron energy, with a 
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FIG. 24: (Color online) The probability for emission of a cer- 
tain number of neutrons as a function of neutron multiplicity 
for thermal neutrons (black circles), En = 4 (red squares), 
9 (green diamonds) and 14 (blue triangles) MeV. The filled 
symbols employ the BSFG level density, while the open sym- 
bols are obtained with a — A/eo. 



total increase of ~ 20% through the energy range shown. 
Although the width of the neutron-energy distribution is 
given by as = E/V2 for a single emission, the resulting 
width grows faster than that with En due to the increased 
occurrence of multiple emissions and thus the appearance 
of spectral components with different degrees of hardness. 

The relatively flat behavior of E{Af) implies that the 
neutron multiplicity ^{Af) will resemble the fragment 
excitation energy E*{Af), as is seen to be the case in the 
bottom panel of Fig.l^where the characteristic sawtooth 
shape of F(A/) is apparent. The number of neutrons 
from the heavy fragment increases somewhat faster with 
En than the number from the light fragment. 

Figure [23] shows the fragment excitation energy E* , 
the kinetic energy of the emitted neutron and the neutron 
multiplicity, all as a function of ^/ for the BSFG and a — 
A/eo at En = 0.5 MeV. The larger x of the BSFG fit gives 
both a stronger dependence of E* on Af and a sharper 
'sawtooth' shape, more consistent with the data shown 
in Fig. [T] The neutron kinetic energy is not strongly 
affected by the value of x. 

New measurements with the fission TPC over a range 
of incident neutron energies could provide a wealth of 
data that could lead to improved modeling. In addition, 
calculations of 'hot' fission that includes temperature- 
dependent shell effects could enhance modeling efforts 
by predicting trends that could be input into FREYA and 
thus test the effects on the PFNS and related quantities. 

For reasons of computational simplicity, we have cho- 
sen to use the same value of x over the entire energy 
range considered. There are some limited data on ther- 
mal neutron-induced fission of ^ssy [32*1 and spontaneous 
fission of ^^^Cf |33|] that suggest the light fragment emits 
more neutrons than the heavy fragment, 40% more for 

U m and 20% more for 252cf [3^. Our BSFG re- 
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FIG. 25: (Color online) The combined kinetic energy of the 
two product nuclei (top) and their residual excitation prior 
to photon emission (bottom) as functions of the neutron mul- 
tiplicity 1/ for thermal neutrons and En ~ 4, 9, 14 MeV. The 
symbols are at the mean values and the vertical bars show 
the dispersions of the respective distributions. 



suit, X ^ 1.23, is consistent with these results. However, 
'hotter' fission could equilibrate the excitation energies of 
the light and heavy fragments which may result in more 
neutron emission from the heavy fragment, also reducing 
the sharpness of the sawtooth pattern. 

Figure |24| shows the neutron multiplicity distribution 
P(i/) for the selected values of En- As expected, V in- 
creases with En and the distribution broadens. However, 
each neutron reduces the excitation energy in the residue 
by not only its kinetic energy (recall E — 2T) but also by 
the separation energy Sn (which is generally significantly 
larger). Therefore the resulting P(i^) is narrower than 
a Poisson distribution with the same average multiplic- 
ity. These results are essentially independent of the level 
density calculation. 

The combined kinetic energy of the two resulting (post- 
evaporation) product nuclei is shown as a function of the 
neutron multiplicity in the top panel of Fig. 1251 It de- 
creases with increasing multiplicity, as one might expect 
on the grounds that the emission of more neutrons tends 
to require more initial excitation energy, thus leaving less 
available for fragment kinetic energy. 
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FIG. 26: (Color online) The combined kinetic energy of the 
two product nuclei (top) and their residual excitation prior 
to photon emission (bottom) as functions of the neutron mul- 
tiplicity u for E„ — 0.5 MeV with the BSFG (squares) and 
a = A/eo (circles). 



prompt neutron emission are not strongly dependent on 
the temperature. 



IV. CONCLUSION 

We have included both multichance fission and pre- 
equilibrium emission into FREYA [1, 0], a Monte-Carlo 
model that simulates fission on an event-by-event basis. 
This has enabled us to perform an extended evaluation 
of the prompt fission neutron spectrum from ^'^^Pu(n, /) 
up to En — 20 MeV. Several physics-motivated model 
parameters have been fitted to the ENDF-B/VII.O eval- 
uation of V and the associated covariance matrix in two 
alternate scenarios for the level-density parameterization. 

Our testing of these two alternate evaluations was in- 
conclusive: neither evaluation performed as well as the 
ENDF/B-VILO evaluation in critical assembly bench- 
marks and, in fact, our two evaluations bracket the 
ENDF-B/VILO evaluation. However, we found improved 
agreement with the LLNL pulsed sphere tests, especially 
just below the 14 MeV peak in the neutron leakage spec- 
trum, for both variants of our evaluation. Although these 
mixed results may limit the utility of our evaluation in 
applications, they do give us hope that further improve- 
ments to the evaluation will either tighten up agreement 
with the critical assemblies or point to other deficiencies 
in the ENDF-B/VII.O 239pu evaluation. 

Further investigations will require fitting to other data 
less sensitive to the v data employed in this work includ- 
ing the albeit low quality PENS data and i>{A) data. 
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The bottom panel of Fig. [25] shows the mass depen- 
dence of the average residual excitation energy in those 
post-evaporation product nuclei. Because energy is avail- 
able for the subsequent photon emission, one may expect 
that the resulting photon multiplicity would display a 
qualitatively similar behavior and thus, in particular, be 
anti-correlated with the neutron multiplicity. 

There is little sensitivity to the calculated level density 
in either shown in Fig. |26| for = 0.5 MeV. 

This result shows that the residual energies left over after 
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